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Internal solitary waves (ISWs) in the NE South China Sea (SCS) are tidally generated at the Luzon Strait. 
Their propagation, evolution, and dissipation processes involve numerous issues still poorly understood. 
Here, a novel method of seismic oceanography capable of capturing oceanic finescale structures is used to 
study ISWs in the slope region of the NE SCS. Near-simultaneous observations of two ISWs were acquired 
using seismic and satellite imaging, and water column measurements. The vertical and horizontal length 
scales of the seismic observed ISWs are around 50 m and 1-2 km, respectively. Wave phase speeds 
calculated from seismic observations, satellite images, and water column data are consistent with each other. 
Observed waveforms and vertical velocities also correspond well with those estimated using KdV theory. 
These results suggest that the seismic method, a new option to oceanographers, can be further applied to 
resolve other important issues related to ISWs. 



Large amplitude nonlinear internal solitary waves (ISWs) are ubiquitous features in the ocean (see www. 
internalwaveatlas.com). The Northeastern (NE) South China Sea (SCS, Fig. 1) is a region where large 
amplitude internal waves are frequently observed\ The observed ISWs in NE SCS are generally considered 
to be generated at the Luzon Strait, the primary water passage connecting the SCS and the Pacific Ocean. They are 
caused by the interaction between the strong tidal forcing and the abrupt topography located there^ l These 
internal waves then propagate westward crossing the deep sea basin region, pass the sharp slope area, entering the 
shoaling shelf zone, and ultimately dissipate on the continental shelf The propagation and evolution processes 
and mechanisms of the ISWs as they cross the SCS are rather complicated due to strong interactions between 
currents, eddies, fronts, internal waves, stratification, and topography^'*'^. These large amplitude ISWs induce 
strong currents which have substantial but not well-studied effects on oceanic mixing, biological processes, 
sediment transport, offshore oil engineering, and submarine navigation 

Although their importance is well recognized, studies of ISWs are far from being comprehensive. In-situ 
hydrographic observations of the ISWs are still rare, and many studies rely only on remote sensing and numerical 
modeling^ '^One new technique that shows promise in observing small-scale features over large areas and full 
depth of the ocean is seismic oceanography, which is capable of capturing finescale structures in the ocean with a 
spatial resolution of O (10 m) using standard seismic reflection profiling^. In less than 10 years, the seismic 
reflection method has been used to image numerous interesting oceanic phenomena, such as internal waves 
subsurface eddies^ horizontal vortices^^, and subsurface currents^^. Nevertheless, the spatial features of near- 
surface ISWs remain unexplored by the seismic reflection method due to the large source-receiver distance of the 
seismic setups and the significant noise of source bubble oscillations. 

In order to image the near- surface water column structures in the NE SCS (Fig. 1), an optimized seismic 
observational system was deployed on 27 June 2012 (Methods). Two prominent near-surface ISWs were clearly 
seen on the seismic image in tandem for the first time (Fig. 2). They were also captured in a satellite image 
obtained for the same day. Propagation velocities of the ISWs derived from the seismic, satellite, and on site 
observation data are all consistent. Our results suggest that the seismic reflection profiling, in conjunction with 
conventional hydrographic observations, could be an excellent method to study near- surface ocean processes 
such as ISWs. 
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Figure 1 | Bathymetry of the NE South China Sea. The purple line shows 
the seismic section carried out from 5:30 to 21:50 on 27 June 2012 (local 
time, UTC + 8:00) with heading of 65 degrees from the north. Red dots 
and the adjacent numbers mark distances from the starting point. Black 
dashed curves are surface signatures of internal solitary waves from satellite 
images^. Red lines are the imaged solitary waves in the study. Figure 
produced with the Generic Mapping Tools (GMT). Bathymetry data is 
from http://topex.ucsd.edu. 

Results 

General features of the water column. A clear seismic image of the 
upper ocean thermocline structure between —50 to 300 m in depth 
was obtained along the continental slope in water depths of 400- 
1000 m (Figs. 1 and 2). Strong features in the images and large 
impedance variations estimated from simultaneous XBT profiles 
are in high correspondence (Supplementary Fig. SI). Further 
comparisons with hydrographic measurements/derivations show 
that temperature is the dominant factor of the seismic response 
(Supplementary Fig. S2). Therefore, it is reasonable to view the 
seismic reflection surfaces as being parallel to isotherms in high- 
vertical gradient areas. 

The most pervasive feature of the seismic image are coherent 
reflections from quasi-horizontal linear features of 5-10 km in 



length, which include undulations with a vertical extent of around 
20 m and a horizontal extent of less than 1 km (Fig. 2a). These are 
likely due to a background field of small internal waves These 
internal waves are more continuous in the deep slope region at the 
right side of the image than in the shallow slope region on the left 
side, where the internal waves are intermittently interrupted by 
patches of blanking zones a few km in horizontal extent and 50- 
100 in vertical extent, possibly indicating regions with well-mixed 
water. The measurement from XBTl within such a patch confirms 
these blanking zones coincide with regions where there is a lack of 
strong gradients in the vertical variation of the hydrographic com- 
ponents (Supplementary Figs. Sla and S2a). The interactions among 
eddy currents, topography, internal waves, or other nonlinear pro- 
cesses are likely responsible for these well-mixed patches. 

Internal solitary waves captured by seismic survey and satellite. In 

addition to the field of small internal waves, two isolated regions 
showing internal vertical displacements of 0(50 m) over a hori- 
zontal scale of less than 1 km were also captured during the cruise 
on 27 June 2012 (Fig. 2). These displacements are interpreted as 
being the large leading waves in packets of internal waves. Each 
wave packet is dominated by one large depression followed by 
smaller waves (Fig. 2b,c). The vertical influence of the initial depres- 
sion extends from 50 m to deeper than 300 m (Fig. 2b). These are 
typical morphological features of ISWs as described elsewhere by 
hydrographic observations or model descriptions^'^^"^^. The ampli- 
tudes of the features are smaller than typically observed for ISWs in 
this region^'^'^^'^^. However, our observations occurred during neap 
tides, when ISWs would have smallest amplitude (Supplementary 
Fig. S3). Hereafter, these two wave packets are named as ISWl and 
ISW2, respectively. Neap tides are also a time when low-frequency 
undular bores have been observed in this region^'^° and such a bore 
may be responsible for the plunging layers between 150 m and 
200 m at the left edge of Figure 2b. 

A NASA Terra satellite image (Fig. 3), acquired on 27 June 2012 by 
the Moderate Resolution Imaging Spectroradiometer (MODIS) 
instrument, shows a number of dark bands forming curving wave- 
fronts with a lateral extent of up to 100 km in the general vicinity of 
the seismic track. These appear to come from two packets, whose 
leading edges cross the seismic track at about —3 and 75 km. By 
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Figure 2 | Seismic image of the water column in the NE SCS. (a) Seismic image of the section in Figure 1 (only the portion of horizontal and 
vertical ranges of 15-70 km and 0-300 m, respectively, is shown). Black triangles are the locations of XBT stations. Detailed structures in the two white 
boxes of (a) are expanded in (b) and (c), respectively. Figure produced with Seismic Unix (SU). 
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Figure 3 | Filtered MODIS image. The observed ISWs at 1 1:05 on 27 June 
2012 (local time, UTC + 8:00) in a MODIS image (http://oceancolor.gsfc. 
nasa.gov). Two red squares mark the locations of ISWl (left) and ISW2 
(right) when observed on the seismic image. Figure produced with the 
Generic Mapping Tools (GMT). 

identifj^ing these wave crests with the two large displacement features 
seen in the seismic section, the space-time differences between the 
seismic and satellite observations can be used to estimate mean phase 
velocities for ISWl and ISW2 to be 1.75 m/s and 2.33 m/s, respect- 
ively (Table 1). 

Mutual verification among seismic, satellite, and site observa- 
tions. In addition to using the space-time difference between the 
seismic and satellite observations (named as approach A), phase 
velocities of ISWl and ISW2 are also obtained using two other 
approaches. These include an analysis of the transient velocities 
from pre-stack seismic data (named as approach B) (Supple- 
mentary Figs. S4-S6)^^'^^ and a theoretical speed derived from 
hydrographic data based on the Korteweg-de Vries (KdV) 
equation (named as approach C)^'^^'^^. 

In approach B, the ISW trough is easy to identify in signals from 
different source -receiver pairs and wave phase speed can be inferred 
from these observations (Supplementary Fig. S4). By assuming a 
constant propagation speed, and aligning the common offset gathers, 
an in-plane propagation speed can be estimated (Supplementary Fig. 
S5). Apparent velocities along the seismic line are 1.88 m/s and 
2.67 m/s for ISWl and ISW2, respectively (Table 1). Then using 
the angle between the ISW propagation and the seismic line esti- 
mated from the satellite image (Table 1; Fig. 3), the apparent velo- 
cities in the seismic line direction are projected into the wave 
propagation direction and the actual wave phase velocities are esti- 
mated to be 1.70 ± 0.06 m/s and 2.19 ± 0.34 m/s for ISWl and 
ISW2, respectively (Table 1). These are consistent with the estimated 
wave speeds using approach A. 

In approach C, the data from XBTl and XBT2 are used to derive 
the analytical parameters of ISWl and ISW2, respectively, because 



the distance of each XBT-ISW pair is very close. With no background 
current and a continuously stratified ocean, the mode-one linear 
wave phase speeds for ISWl and ISW2 are estimated as 1.40 and 
1.68 m/s, respectively (Table 1). These wave speeds can be corrected 
for nonlinear effects by assuming that a continuously stratified 
Korteweg-de Vries equation adequately describes the propagation. 
The nonlinear phase speeds are 1.63 m/s and 1.92 m/s. Here, the 
estimated phase speed of ISWl is comparable to the observed values 
while the speed of ISW2 is about 12% lower than the seismic obser- 
vation, although still within the uncertainty of the seismic obser- 
vation. This discrepancy between the theoretical estimates and the 
observed phase speeds may be caused by the neglected background 
currents^^'^^ although satellite altimetry- derived velocity fields (not 
shown) suggest that the mean current at least is approximately per- 
pendicular to the direction of propagation. Another explanation for 
the difference is that ISW2 appears to be transforming its shape and 
is not a pure, symmetric solitary wave. 

Other important parameters during the wave propagation stage 
such as waveforms and vertical velocities can also be derived from the 
seismic data and compared with analytical estimates from KdV 
theory. 

While the seismic image directly provides the waveform of ISW by 
tracking the reflective horizons, the hydrographic data can also pre- 
dict the wave shape using the theoretical models. In this study, the 
ratios of the amplitude to the upper-layer depth are around 0.7, so a 
weakly nonlinear theory^^, e.g., the KdV theory may be sufficient. 
The analytical solitary solution of the classical KdV equation 

— +c—+(xrj—+P —-r =Ois r](x, t) = rjosech\(x - vt)/A), where 

dt ox ox ox^ 

the nonlinear velocity v and the characteristic width A are related to 
the mode-one linear speed c, the nonlinear parameter a, the disper- 
sion parameter jS, and the maximum wave amplitude rj^hy v = c + 
af]o/3, = 12j5/af/o (Supplementary Table SI). The ratios of the 
seismic observed full widths at the half wave amplitude to the char- 
acteristic widths are 1.71 and 1.73 for ISWl and ISW2, respectively, 
while the analytical ratio is 1.76. Therefore, the KdV theory predicted 
solitary wave shapes match well with the seismic observed wave- 
forms (Fig. 4). 

Similarly, the vertical velocities of the observed ISWs are also 
compared with analytical predictions of the KdV theory. From the 
seismic image, the vertical velocity between the peak and the half 
peak of the displacements is measured using the relation w = rjo/L 
(Supplementary Fig. S7). Estimated vertical velocities are 13.8 cm/s 
and 8.7 cm/s for ISWl and ISW2, respectively. In comparison, KdV 
theory predicted values are 12.9 cm/s and 7.5 cm/s using the equa- 

tion w(x)= — -^^^t=o= —2-^sech^(x)tanh^(x) (Supplementary 
Fig. S8). 

Discussion 

Seismic profiling has previously been used to image oceanographic 
features at depths of 100 s to 1000 s of meters in the ocean. Less 
attention has been paid to features closer to the surface, because the 



Table 1 | Parameters observed /derived during velocity estimation. Tsat: satellite acquisition time; Tseis: seismic acquisition time; LT: local 
time, equals to UTC + 8:00; Dist: ISW distance between satellite and seismic observation; A: angle between ISW propagation and seismic 
line heading; rj^: amplitude of ISW. Vsat: velocity from satellite observation; Vsa: apparent phase velocity from seismic data along the 
seismic line; Vseis: phase velocity from seismic data along the wave propagation direction; c: mode-one linear phase speed; KdV: speed 
estimated using continuous stratified ocean system; Vvs: vertical mean velocity from seismic data; Vvx: vertical mean velocity from 
hydrographic data 

Tsat(LT) Tseis (LT) Dist (km) A (°) ^Jm) Vsat(ms-^) Vsa(ms-i) Vseis (ms"^) c{ms-') KdV(ms-^) Vvs (cm s'^) Vvx (cm s"^) 

ISWl 11:05 07:45 21.0 25 60 1.75 1.88 1.70 1.40 1.63 13.8 12.9 

ISW2 11:05 12:19 10.5 35 50 2.33 2.67 2.19 1.68 1.92 8.7 7.5 
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Figure 4 | Waveform comparison of ISWs. Two predicted waveforms of ISWs (white lines, derived using the stratified KdV theory from XBTl and XBT2 
profiles) are superimposed on their seismic images. The predicted waveforms have been contracted by simulating the Doppler-Effect for wave-ship 
movements and then projected to the direction of seismic section. Figure produced with Seismic Unix (SU). 



conventional geometries of marine seismic survey are unfavorable 
for imaging regions that close to the surface. However, by optimizing 
the survey geometry, we have successfully formed a seismic image of 
the upper 300 m of the water column (but below the mixed layer) in 
the NE SCS. The image shows that thermocline fine -structure in this 
area is dominated by typical internal waves with a vertical extent of 
0(10 m) that can stretch for up to 10 km. These features are wea- 
kened or interrupted by patches a few km across where vertical 
profiles are smoother, verified by in-situ XBT observations. The most 
prominent features are the two distinct displacements in the section, 
associated with ISWs that form discrete packets radiating away from 
their source regions during neap tide when ISWs are rarely observed 
in this area^'^. The in-plane propagation speeds of these displace- 
ments were estimated from the pre-stack seismic data alone, and 
the results are consistent with propagation speeds estimated in other 
ways. In addition, parameters like wave shapes and vertical velocities 
are successfully determined from the seismic data, and can be com- 
parable to the theoretical predictions. 

Thus we have shown that studying ISWs becomes more compre- 
hensive using an optimized seismic array and incorporating conven- 
tional observations. Other fme-structure features observed on the 
seismic image await confirmation. One of these features are the 
intriguing regions of possibly enhanced mixing occurring immedi- 
ately after passing of the ISWs, identified by clear disappearance/ 
weakening of the reflections around 200 m depth (Fig. 2), which 
may be a result of shear instability^^. Although weak stratification 
in this area is more easily destroyed than the strong stratification 
above, whether such appearances are universal and whether they 
represent quantitative descriptions of mixing is unknown. Another 
feature requiring explanation is the flat-step structure of the ISW2 
trough which is unable to be modeled by single KdV soliton. Its 
dynamic implication (wave splitting or merging, or other process) 
also remains to be assessed. 

Therefore, we believe that continued studies of other interesting 
stages or phenomena of ISWs, such as generation, evolution, diffrac- 
tion, dissipation, polarity conversion, and mixing will benefit from 
coordinating seismic observations with other conventional 
techniques. 

Methods 

Seismic observation and processing. A seismic oceanographic experiment was 
conducted aboard the R/V SHIYAN2 in the NE South China Sea on 27 June 2012. In 
order to optimize for imaging near- surface features, a Sercel GI gun was chosen as the 
seismic source since it has good performance in suppressing bubble oscillations. 
However, due to the limited source energy from a single Sercel GI gun (with volume of 
105/105 cu-inch), a modified acquisition geometry was necessary to improve the data 
quality. The nearest offset (source-receiver distance) was set to 90 m to reveal the 
shallow structures, which are overlooked by the conventional geometry in which the 



nearest offset is around 250 m. A Sercel Seal oil filled streamer with 72 receivers (12.5- 
m interval) was towed at 5 m depth. The air gun was triggered every 25 m (~ 10 s) to 
increase the fold number. The record length was 6 s with sampling interval of 2 ms. 
Such a geometry ensured a satisfactory ray sampling of the near- surface structure and 
a high SNR (Signal Noise Ratio) of the stacked seismic image in processing. 

Seismic data were processed using the pre-stack depth migration (PSDM) 
method^*'^^ The conventional procedure of PSDM consists of six steps: (1) geometry 
definition, (2) bad trace deletion, (3) direct wave attenuation, (4) band-pass filtering, 
(5) PSDM with velocity analysis, and (6) common image gather stacking. Here, the 
velocity analysis of the fifth step was omitted. Instead, the acoustic velocity was fixed 
at 1510 m/s according to our hydrographic observations. Migration tests have shown 
that it is an acceptable model with fine imaging. According to the seismic geometry 
and processing parameters in the study, the shallowest reliable seismic reflection is 
—50 m depth, which is similar to the mixed layer depth from the XBT data. 

Velocity estimation of the ISWs from seismic data. The propagation velocities of 
the ISWs can be estimated using the pre-stack seismic data because of the redundancy 
of the seismic acquisition spanning a finite time^^'^^'^*^. However, the near-surface 
reflection signals from ISWs tend to be muted by normal move out (NMO) and 
suppressed by direct wave attenuation. Thus the size of the common mid-point 
gathers (CMPs) or folds of the datasets decreased from 18 to around 7 traces (i.e., only 
the first 27 near-source traces of the common offset gathers (COGs) are left over). 
Therefore, the time and distance spans are too short to estimate the movement of the 
reflections precisely and intuitively according to Sheen's approach using the CMPs^'*. 
To overcome this limitation, we resort to the remained COGs which can be used to 
track the ISWs clearly and improve the approach accordingly. Integrity of the 
waveforms on the COGs helps us to derive the crest/trough points used for movement 
estimation. 

The principle of the in-plane (apparent) velocity measurement is illustrated in 
Supplementary Fig. S4 with the assumption of uniform vessel speed. For the given 
coordinate system, horizontal locations are, cmp2 = — (g2*dg + ofs0)/2, cmpl = 
— (gl*dg + ofs0)/2 — (s2 — sl)*ds, where(sl,gl) and(sl,g2) are the numbers of the 
shot- receiver pairs that capture the crest/ trough of the moving internal wave at the 
midpoints cmpl and cmp2, respectively; dg, ds, ofsO are the receiver interval, shot 
interval, and nearest offset, respectively. Thus, the horizontal velocity of the wave 
crest is measured by v = (cmp2 — cmpl)/T = — [(g2 — gl)*dg/2 — (s2 — sl)*ds]/ 
[(s2 — si )*dt], where T is the shot time duration between si and s2, and dt is the shot 
time interval. Taking ISWl as an example, we first auto-tracked the waveforms from 
27 COGs, and then interpolated to find the troughs of the ISW and the corresponding 
source-receiver pairs. Then, with the known receiver span (g2 — gl) of 26, the only 
unknown shot span (s2 — si) was estimated using the slope of the regression based on 
the RMS criterion (Supplementary Fig. S5) . Finally, the apparent horizontal velocity v 
was derived, as well as the actual phase velocity on the direction of wave propagation. 
For the real data, the source range, time duration, and ship heading were estimated 
precisely using the navigation data of the cruise. 

To confirm the utility of this velocity estimation method, we estimated the velocity 
of a static seamount, whose spatial scale was comparable to ISWs, to perform the 
whole procedure to measure its velocity (Supplementary Fig. S6). The estimated 
velocity for the static seamount was only -0.04 m/s, suggesting any bias is small. A 
larger source of uncertainty is from the slope regressions and 95% confidence 
intervals estimated from the regression were ±0.06 m/s and ±0.34 m/s for ISWl and 
ISW2, respectively. The primary factor of the large error of ISW2 is that the troughs of 
the ISW2 on COGs are too gentle to allow the trough to be accurately located. 

Velocity estimation of the ISWs using XBT data. We used the in situ observations 
from XBTl and XBT2 to calculate the nonlinear phase velocities of ISWl and ISW2 
correspondingly. Here, the salinities of the profiles were estimated using the historical 
CTD data in the study region from the World Ocean Database 2009^^. Thus, the 
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density, velocity, impedance, buoyancy frequency were derived using the 
temperature/salinity/depth relationships. For a continuously stratified fluid, the first- 
mode linear velocities c were solved from the Sturm-Liouville equation^^ 
(fWiz) 

--\-^ W(z) = 0, W(0) = W(H) = 0. Then, the parameters a and of the one- 



drj 



-- 0 were calculated for the 



dimensional KdV equation 

same stratification^^, as well as the characteristic width A (Supplementary Table Sl)^. 
Finally, the nonlinear phase velocities were estimated via v = c + (xrjQ/3. 



dri dri ^d^ri 
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